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Non-relativistic collisionless shock waves are widespread in space and astrophysical plasmas and are known 
as efficient particle accelerators. However, our understanding of collisionless shocks, including their structure 
and the mechanisms whereby they accelerate particles remains incomplete. We present here the results of 
numerical modeling of an ion-acoustic collisionless shock based on one-dimensional (ID) kinetic approximation 
both for electrons and ions with a real mass ratio. Special emphasis is made on the shock-reflected ions as 
the main driver of shock dissipation. The reflection efficiency, velocity distribution of reflected particles 
and the shock electrostatic structure are studied in terms of the shock parameters. Applications to particle 
acceleration in geophysical and astrophysical shocks are discussed. 
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I. INTRODUCTION 

A fundamental problem of dissipation of flows in colli¬ 
sionless plasmas moving faster than all plasma modes 
available has numerous laboratory, astrophysical and 
geophysical applications. It was realized more than half a 
century ago (see e.g.^^^^^ for reviews) that despite lacking 
collisions, a shock wave must form in such flows. The new 
applications pose now new challenges to our understand¬ 
ing of collisionless shocks structure and their potential to 
accelerate particle^^^^. Microscopically, strong shocks 
can be supported by suitable plasma waves that ran¬ 
domize particle trajectories in lieu of binary collisions. 
Although these waves are not necessarily energetically 
important, in many cases they are driven by shock re¬ 
flected particles. In this paper we consider relatively 
weak shocks in which the crucial process of ion reflection 
can be effectively approached from the well understood 
physics of solitary waves. 

The shock at first place results from a nonlinear steep¬ 
ening of a certain mode whose dispersion often limits the 
steepening. For a number of modes nonlinearity and dis¬ 
persion can balance each other exactly which results in 
a soliton. In this paper we consider the most fundamen¬ 
tal ion-acoustic solitons with a trailing wave trairP^^, 
which are the building blocks of collisionless shock waves 
in nonisothermal plasmas, with the electron temperature 
much higher than the ion temperature, Tg ^ T^. Our 
study is relevant to the elect rostatic shock propagation 
in laser-produced plasma^^ElHIIl, especially to the gener¬ 
ation of monoenergetic ion beams, and to shock-related 
processes in astrophysical and space plasmas. Among 
astrophysical applications of this study we would like to 
mention the ion injection into the first order Fermi accel¬ 
eration in astrophysical shocks. The injection process is 
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crucial for the resolution of such problems as the prob¬ 
lem of cosmic ray origin. Overall, the cosmic ray spec¬ 
trum is a power-law starting from supra-thermal energies 
and the injection efficiency controls the normalization of 
the spectrum. It ensures a smooth transition between 
the thermal plasma and the power-law spectrum at high 
energies. Although the power-law component is due to 
multiple shock front crossing and not so much due to the 
reflection from the shock, the particle reflection can also 
contribute to the high energy part of the spectrum, if the 
shock is oblique. In this case the magnetic field is impor¬ 
tant. The electrostatic mechanism is not efficient at high 
particle energies. However, our purely electrostatic simu¬ 
lations are relevant to the microphysics of particle reflec¬ 
tion off the shock front including the cosmic ray loaded 
shocks. 


II. MODEL AND SIMULATION SET-UP 

In order to address the ions reflected and escaped from 
the shock, incorporate the ion reflection into the global 
shock structure and investigate its effect on the shock it¬ 
self we use the kinetic plasma description. The plasma 
dynamics is governed by the Vlasov equations for the dis¬ 
tribution functions of the plasma components and Pois¬ 
son equation for the electrostatic field. Numerically, we 
solve the Vlasov equations using Particle-in-cell (PIC) 
methoeP^, the electrostatic field is obtained by direct in¬ 
tegration of the charge density, which in turn is calculated 
from the individual particle positions. For the solution 
of the equations of motion of particles an organic second 
order in time and space time-reversible numerical algo¬ 
rithm is applied. 

Fig. g shows the schematic set-up for our ID simula¬ 
tions. A shock wave is produced by reflection of a high¬ 
speed electron-ion plasma off a conducting wall. The 
shock forms due to the interaction between the incom¬ 
ing and reflected flows and propagates away from the 
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FIG. 1. Scheme of the set-up for ID kinetic modeling of elec¬ 
trostatic shock wave 

wall. We measure the density of plasma components 
in units of the unperturbed density ng, coordinates in 
units of the Debye length of the unperturbed plasma flow 
Xd = \/Te/ 47 re^no, velocity in units of the electron ther¬ 
mal speed uq = y/TXjnie, and time and electric field in 
units of l/wpe = -yTOe/47re^no) and Eq = V 47 rnoTe cor¬ 
respondingly. 

In out PIC simulations, we ensure the high enough 
resolution to resolve low-density regions, sharp gradients, 
and the dynamics of both electrons and ions by taking at 
least 100 particles per cell and a spatial resolution bet¬ 
ter than 0.05 Xd- For reference, in all simulations shown 
the conducting wall is located at a: = 0 and the plasma 
flow impinges on it at t = 0. Within these settings, the 
dimensionless input parameter which determines the evo¬ 
lution of the plasma flows is Mg = Vb/Cs, where Vq is 
the speed of unperturbed plasma flow and Cs = i/Te/rui 
is the sound speed. Moreover, the stability of the soliton 
and the formation of the shock are proven to be strongly 
dependent on the velocity distribution of ionJi^, thus on 
the fraction Ti/Tg. 

Note, that this purely electrostatic approach is reason¬ 
able for quasi-parallel shocks as it treats scales Xd < 
I < cjwpi ~ r* (as ,0 ~ 1 very often), thus providing 
an important microphysics ingredient for the diffusive 
shock acceleration as welP®^. The hybrid simulations, 
broadly used in astro-geo-shocks studies, do not resolve 
these scales since they treat electrons as a fluid. 

III. SIMULATION RESULTS 
A. Cold upstream ions 

The evolution of the distribution of ions on phase plane 
(x, 14) for Mg = Vg/Cs = 1.3 is shown in Fig. a)- 
d). By t = 200 ^ 4 ,^ the flows interacted and the shock is 
formed. The isolated density (and electrostatic potential) 
structures. Fig. e)-f), have the characteristic width 
of several Xd- Ions reflected from the shock are clearly 
visible on the phase plane. In the case of cold upstream 
ions the ions are reflected at a single point where the 


potential reaches its maximum and the electric field has 
a jump and the potential profile a cusp. 

The number of ions reflected out of the shock increases 
with the increasing of the flow velocity, see Fig. where 
the distributions of ions on the phase plane (X, I 4 ) and 
the distribution functions of ions f{v) at tujpe = 2800 
are shown. The distribution function has three distinct 
peaks: inflowing ions {v = Vg), reflected ions (|z)| ~ (Vg-f 

front), with Vfront the Speed of the reflecting hump), 
and ions at rest downstream of the shock. The energy 
spectrum of the reflected ions is very broad. 

The spatial distributions of the electrostatic potential 
{—ecj)) for different values of upstream flow are presented 
in Fig. [^a)-b). It is clearly seen that when the ions begin 
to reflect from the soliton tip, the classical single solution 
bifurcates into a more complex structure that comprises 
the leading soliton, the periodic wave train downstream 
of it and the foot occupied by the reflected ions. This foot 
is supported by the reflected ions and also accelerates 
them somewhat further. Because of the reflection of cold 
ions the leading soliton for Mg = 1.3, 2.1 and Mg = 3 
acquire a clear cusp structure. 

The particle reflection strongly depends of the ampli¬ 
tude of the shock (overshoot of electrostatic potential). 
In turn it alters the shock amplitude and speed, and thus 
the reflection threshold. Moreover, the reflected ions per¬ 
turb the electrostatic potential upstream which changes 
the speed of inflowing ions and thus, again, the condition 
for the ensuing reflection. For Mg = 2.1 and Mg = 3 
the shock looses much of its momentum because of the 
ion reflection, for Mg = 1.3 the number of reflected ions 
is very small. Momentum conservation suggests that ion 
reflection slows down the shock, thus decreasing its am¬ 
plitude. Evidently, decreasing amplitude acts against re¬ 
flection. However, if the amplitude drops below the re¬ 
flection threshold for a given speed, the shock amplitude 
should recover and the reflection will resume. By virtue 
of these feedback loops the reflection efficiency can be 
expected to either stay at a self-organized critical level 
or to exhibit cyclic dynamics. Numerical simulations and 
observations validate both scenarios, though, the cyclic 
shock behavior was prim arily observed in magnetized 
plasma environments^^BlI. The temporal evolution of the 
maximum value {—ofmax) of the electrostatic potential 
extracted from our simulations. Fig. c) (solid line), 
confirms a cyclic behavior of the shock. 

The approximate velocity of the shock with respect 
to the upstream flow can be defined from the position 
of the maximum of the electrostatic potential (dashed 
curves in Fig. [^c)). The result is shown in Fig. Eld). 
Measured in units of Cg, this velocity directly represents 
the shock Mach number Ma = Vshock/Cs- Our simula¬ 
tions prove that the solution, consisting of the leading 
soliton, periodic wave train downstream and the foot oc¬ 
cupied by the reflected ions persists with the increasing 
Mach numbers up to a critical value of Ma = M 2 ~ 4.5 
until almost all incident ions get reflected. In case of cold 
ions this critical value is reached for the upstream veloc- 
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FIG. 2. a)-d) - The distribution of the ions on the {X, 14) phase plane. e)-f) - The distributions of the electrostatic potential 
and ion density att = 880,1520, 2160, 2800 for Mq = VofCs = 1.3. The leading soliton has a cusp structure. It occurs only 
when the cold ions get reflected. 



FIG. 3. a)-d) - Distribution of ions on the (X,Vx) phase plane for different values of Mq = Vo/Cg. e) - Distribution function 
of ions fi{vi). 


ity Vb = V 2 = S.dSCs. The obtained critical value of M 2 
is much beyond a critical value M*, at which the shock 
is about to reflect some of the upstream ions. For the 
cold upstream ions and Bolzmannian electrons (this is 
definitely not the case in our simulations) M* ~ l.lP, for 
adiabatically trapped electrons M* ~ 3.1^. Note, that 
in the early days of collisionless shock research the exis¬ 
tence of stationary shock transitions for Mach numbers 
exceeding critical M* has been often disfavored. 

When the velocity of the incoming flow is bigger than 
V 2 the solution becomes qualitatively very different. An 
almost monoenergetic beam of ions, produced by the re¬ 
flection of the inflowing ions with Mq = 3.85, is clearly 


seen in Figj^ Note also the sharp peak in the distri¬ 
bution function (vi < 0) and the significant deceleration 
of the inflowing ions {vi > 0) for Mq = 3.85. The next 
qualitative change in the solution should occur when the 
inflow velocity reaches the values comparable to the ther¬ 
mal velocity of the electrons. At this point the Buneman 
instability starts to play an important role in the interac¬ 
tion of flows. This regime is however outside the current 
study. 

Finally, frame a) of Fig. shows the distribution of 
ion density ni{x, t) in logarithmic scale and the time his¬ 
tory Xi{t) of some selected ions. The first peak in the ion 
density corresponds to the “outrunner” soliton, liberated 
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FIG. 4. a)-b) - Electrostatic potential for Mq = 0.4,1.3, 2.1 and 3. c)-d) - Evolution of the maximum of the electrostatic 
potential and the velocity of the shock wave in time. The cusp structure in the potential occurs if cold ions are reflected. The 
outer soliton is liberated from its dispersive tail. The solution is asymmetric with respect to the reflection point, its downstream 
is oscillatory. Upstream of the soliton, reflected ions create a foot with an elevated electrostatic potential. 




FIG. 5. (color online) Distribution of ions on (X,Vx) phase 
plane (left frames) and distribution function of ions fi(vi) 
(right frames) at tujpe = 2800 for Mo = 3.43 and Mq = 3.85. 
For Mq = 3.85 the flow velocity of the incident plasma exceeds 
the critical value V 2 . 


from its dispersive trail. The velocity of this soliton is al¬ 
most constant in time (apart from the oscillations, men¬ 
tioned before) and it reflects a lot of incoming ions, which 
then escape with double the shock speed, see also Fig. 
b). The trail, instead, not only propagates with a smaller 


velocity, it slows down significantly as well. Some of the 
ions get trapped between the soliton and the leading part 
of the trail. 


B. Warm upstream ions 

The resolution of our PIC simulations is high enough 
to accurately include the finite ion temperature in our 
analysis. Fig[^ shows the distributions of ions on the 
phase plane for Mq = 0.4 (very slow upstream flow) and 
cold/warm upstream ions. If the upstream ions are cold 
{Ti = 0) only a negligibly small number of them get re¬ 
flected, but if Ti > 0 the final distribution function of 
ions assumes a different shape. In Figj^ a steady reflec¬ 
tion of ions is clearly visible for T^ > 0. Their number 
remains, however, very small and the energy spectrum 
- quite broad. Moreover, there are ions trapped in the 
“multipeaks”. Note, that stability of solitons and forma¬ 
tion of the electrostatic shock waves in ultraintense laser 
interaction with overdense plasmas are also found to be 
strongly dependent on the velocity distribution of ionJi^. 
It should be mentioned, that if the flow velocity of warm 
(both electrons and ions) inflowing plasma exceeds V 2 the 
solution qualitatively similar to Figj^forms again: the re¬ 
flection produces a huge amount of almost monoenergetic 
ions which propagate away from the shock. 
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FIG. 6. a) - Ion density distribution and ion time history for Mq = 2.14 b) - Particles trajectories on {X{t), I4(t)) plane. 



FIG. 7. Distribution of the ions on {X, 14) phase plane for Mq = Vo/Cs = 0.4 and Ti = 0 (left frame), Ti = 0.045re (middle 
frame), Ti = O.OOTe (right frame). 


IV. CONCLUSION 

Our simulations show that (i) when the Mach num¬ 
ber of the acoustic soliton increases to the point of ion 
reflection and the soliton transforms into a soliton train 
downstream, this structure persists with increasing Mach 
number until almost all incident ions get reflected; (ii) the 
ion reflection alters the shock amplitude and speed, thus 
impacting the reflection threshold itself; (iii) the pedestal 
of the electrostatic potential, supported by the reflected 
ions, changes the speed of inflowing ions and thus again, 
the condition for the ensuing reflection. The formation 
of the pedestal in the potential between the reflecting 
hump and the front edge of the group of reflected ions 
makes the analytic theory very challenging. In fact, this 
pedestal is sensitive to the model used for the transi¬ 
tion between the reflected ions and upstream (e.g. cold 
vs hot beam of reflected ions). In the pedestal the in¬ 
flowing plasma is slowed down considerably, if its front 
edge moves at the speed of the main shock, but reflected 
particles pass through it and accelerate, decreasing their 
density. All these uncertainties are critical and hard to 
resolve without simulations. Note, that an analytic so¬ 
lution for the ion-acoustic collisionless shock with self- 
consistently reflected ions and their subsequent cooling 


was recently presentecP^^^. Its relation to the well known 
non-reflecting ion-acoustic soliton solution, strictly lim¬ 
ited by the critical Mach number M* = 1.6 (Boltzmann 
electrons )P and M* = 3.1 (trapped electrons)!^, was clar¬ 
ified and quantified. 
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